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Abstract 

We study associative memory neural networks of the Hodgkin-Huxley type of spiking 
neurons in which multiple periodic spatio-temporal patterns of spike timing are memo- 
rized as limit-cycle-type attractors. In encoding the spatio-temporal patterns, we assume 
the spike-timing-dependent synaptic plasticity with the asymmetric time window. Anal- 
ysis for periodic solution of retrieval state reveals that if the area of the negative part of 
the time window is equivalent to the positive part, then crosstalk among encoded patterns 
vanishes. Phase transition due to the loss of the stability of periodic solution is observed 
when we assume fast a-function for direct interaction among neurons. In order to eval- 
uate the critical point of this phase transition, we employ Floquet theory in which the 
stability problem of the infinite number of spiking neurons interacting with a-function is 
reduced into the eigenvalue problem with the finite size of matrix. Numerical integration 
of the single-body dynamics yields the explicit value of the matrix, which enables us to 
determine the critical point of the phase transition with a high degree of precision. 

1 Introduction 

Synchronized firing of neurons are ubiquitous phenomena in real nervous system, and capa- 
bility of synchronicity of neurons for information processing has been the subject of many 
research papers . It has been revealed that repeating firing patterns of 

pyramidal neurons appear in sharp waves of rat hippocampus [[l^]. This result of experiment 
suggests the possible role of spatio-temporal patterns of spike timing in encoding information 
in real nervous system. Associative memory neural networks that memorize spatio-temporal 
patterns of spike timing is essential for understanding this information processing of spike 
timing. 

Much of fundamental concepts of associative memory neural networks have been devel- 
oped by replica calculation of Ising spin neural networks with the energy function [[I3l|l4[|l5|]. 
In these neural networks, the standard type of Hebb rule is assumed to define symmetric 
synaptic connections, which bring about fixed-point-type attractors. These fixed-point-type 
attractors are, however, useless for encoding spatio-temporal patterns. Asymmetric synaptic 
connections play a significant role in encoding spatio-temporal patterns, and then the ques- 
tion arises about the learning rule that defines asymmetric synaptic connections so that the 
network functions as associative memory for sptio-temporal patterns. When we assume syn- 
chronous update rule for the dynamics of spin neural networks, a simple extension of the 



Hebb rule readily realizes associative memory for spaio-temporal patterns |16|. Neverthe 



less, the problem becomes rather difficult when we assume asynchrnous update rule for spin 
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neural networks. Complicated learning rules are required to control the continuous transition 
of network state in sequential retrieval of spatial patterns [jT^,[T8|]. 

In spin neural networks | ]l9| , p0| , pl] ], as well as analog neural networks [ p2| , p3| , p4| , p5| , p6| , p7| 
^^9ll, state variables of neurons are assumed to represent their firing rate. In neural networks 
of phase oscillators, phase variables are used to represent synchrnoized firing of neurons. 
The synaptic connections of Hermitian permits networks of oscillators to memorize spatio- 
temporal patterns of phase differences. Since some theoretical techniques are available for the 
analysis of phase oscillators, the properties of networks of oscillators have been investigated 
extensively [ ^pl] , ^^|34l , ^ ] . Even in the presence of white noise as well as heterogeneity 
of oscillators we can derive the storage capacity of networks of oscillators analytically [^]. 

Nevertheless, networks of oscillators may possibly offer a distorted interpretation of syn- 
chronized firing in the real nervous system unless interactions among neurons are sufficiently 
weak. To provide a real understanding of the information processing of spike timing, we 
must adopt more biologically plausible models of neural networks. For this purpose, neural 
networks of spiking neurons are considered to be suitable models for investigation, though it 
remains an unsolved problem to find adequate learning rule for spatio-temporal patterns of 
spike timing. Since asymmetric synaptic connections bring about sequential firings of spik- 
ing neurons ||3^,|37|], one may consider that asymmetric synaptic connetions are essential for 
associative memory neural networks of spiking neurons. In fact, incorporating asymmetric 
synaptic connections, Gerstner et al. have successed in encoding a few number of spatio- 
temporal patterns in networks of spiking neurons with discrete time dynamics [^]. 

The spike-timing-dependent synaptic plasticity found in electrophysiological experiments 
excites a good deal of interest in this connection. It has been revealed that the modification of 
excitatory synaptic weight depends on the precise timings of presynaptic and postsynaptic fir- 



ings [|38|,p9, 40|. Synaptic weight is found to increase if firing of a presynaptic neuron occurs 
in advance of firing of a postsynaptic neuron, and to decrease otherwise. The spike-timing- 
dependent synaptic plasticity is described by the time window with the negative part as well 
as the positive part (Fig. [l|), and this asymmetric shape of the time window has been attract- 
ing a growing interest of reseachers [|4l|,|4^,p3|,p4],^ p^|4^ ,^]. Since the asymmetric time 
window brings about asymmetric synaptic connections, the spike-timing-dependent synaptic 
plasticity is thought to be advantageous to encode spatio-temporal patterns. In the previous 
study of ours, we have studied associative memory neural networks of spiking neurons in 
which the asymmetric time window of the spike-timing-dependent synaptic plasticity is used 
to encode multiple periodic spatio-temporal patterns of spike timing [[lO[]. We have assumed 
networks of the Hodgkin-Huxley neurons interacting through direct synaptic interaction, as 
well as indirect synaptic interaction intermediated by firings of interneurons. In the process 
of memory retrieval, the indirect interactions bring about the oscillatory inhibitory electric 
currents, which regulate spike timings of neurons as in the case of gamma and ripple oscil- 
lations ^ ^1]] . In order to elucidate the stationary properties of these retrieval state we 
have derived the periodic solution for retrieval state analytically, and then we have shown that 
if the area of the negative part of the time window is equivalent to the positive part, crosstalk 
among encoded patterns vanishes. This theoretical result implies the outstanding nature of the 
spike-timing-dependent synaptic plasticity for encoding multiple spatio-temporal patterns. 

In our previous study, however, we did not carry out a stability analysis for the retrieval 
state, and hence it remained unclear whether the derived retrieval state are stable or not. We 
investigate the same models of neural networks also in the present study, but we assume that 
the a-function of the direct interaction decays much faster than the previous one. Then, we 
find the phase transition due to the loss of the stability of retrieval state. In order to evaluate 
the critical point of this phase transition we consider employing Floquet thoery. Neverth- 
less, the degree of freedom of the present system is infinite, and the naive application of 
Floquet theory yields the eigenvalue problem with the infinite size of matrix. Furtheremore, 
a-function we assume here exhibits the infinite long-time influence, and its treatment may 
also require the infinite size of matrix [|2|l. Without calculating the explicit form of the ma- 
trix, Bressloff et al. have investigated the stability of some periodic solutions in networks 
of integrate-and-fire neurons [|3]|, though its application to other neuron models, such as the 
FitzHugh-Nagumo neurons and the Hodgkin-Huxley neurons, seems to be limited. 
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In the present study, we employ two theoretical techniques to reduce the size of the matrix 
for Floquet theory. We first take the limit of the infinite large number of neurons, which 
reduces the stability problem of neurons into the problem of Q sublattices. Then, we define 
some additional variables for each sublattice and evalute the infinite long-time influence of 
a-function with the finite size of the matrix. The explicit value of the matrix is calculated 
from the numerical integration of the single-body dynamics. Therefore, we can explicitly 
obtain the eigenvalues of the matrix, which enable us to determine the critical point of the 
phase transition with a high degree of precision, even when we assume networks of Hodgkin- 
Huxley neurons. 

The present paper is organized as follows. In section ^ we present the dynamics of 
neural networks of spiking neurons, and then introduce the spike-timing-dependent learning 
rule for associative memory. In section ||, we derive the retrieval state analytically in the limit 
of the infinite number of neurons. After that, the stability of this retrieval state is analyzed 
by Floquet theory in section ^ In section ^, we illustrate the typical behavior of network 
in the process of memory retrieval. The result of numerical simulations are presented, and 
compared with the theoretical results. In section |[ we discuss the phase transition due to the 
loss of the stability based on the stability analysis in section Q In section ^ we investigate the 
case of slow a-function, with which we find two separated retrieval phases. In one of these 
retrieval phases, neurons obtain the large size of the oscillatory inhibitory synaptic electric 
currents, which well regulate the spike timing of neurons. Finally, in section ||, we give a 
brief summary and discuss the biological implication of the present study. 



2 Associative memory neural networks of spiking neurons 
2.1 Network dynamics 

In real nervous system, many regions such as the neocortex and hippocampus are found to 
comprise a large number of pyramidal neurons as well as interneurons. Our interest in the 
present study lies in spike timing of pyramidal neurons, and we denote the dynamics of 
pyramidal neurons by a set of nonlinear differentail equations of the form 

V/ = /(v,.,w,.i,...,w,.„) +/,-, (1) 
^,7 = ^/(^P%p---'^m)> / = !,...,«,/ = 1,...,A^ (2) 

with 

A' = ^PP,; + V + ^EXT.i ' 

where v,- denotes membrane potential and auxiliary variables w.^ are used to describe gating 
of ion channels. Synaptic electric currents denote interaction among neurons, and the defi- 
nitions of three currents included in will be given later For the dynamics /(v^jW^j , . . . jW^^) 
and (I'p , . . . , w,.), many authors assume the integrate and fire equation, the FitzHugh- 
Nagumo equations [p4|,|55[], the Hodgkin-Huxley equations |]56|], and so on. In the present 
study we choose the Hodgkin-Huxley equations, which are summarized in Appendix^ 

The synaptic electric current /pp ■ in expresses the direct interaction among pyramidal 
neurons. We define spike timing of neuron / as the time when the membrane potential 
exceeds the threshold value = 0. Denoting ^-th spike timing of neuron / by r, (A:), we define 
the current /pp ^ as 

7=1 k 

where J^^ represents the synaptic weight, and a-function 5pp(f ) is defined as 

To for f < 0, 

'^pp(0 = { \ /'g-'/Tppi _e-'/-rpp.2'j for < f. 

\ '^VV.X ^ '^PP,2 ^ ' 
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The constant App is used to control the intensity of the current Ipp ■. In the following section, 
we will investigate the case of the fast a-function 5'pp(r) with Tpp j = 3 [msec] and Tpp2 = 
0.3 [msec] as well as the slow a-function Spp{t) with Tpp j = 20 [msec] and Tpp j = 2 [msec]. 

The synaptic electric current /jp in /■ expresses the indirect interaction among pyramidal 
neurons that is intermediated by firings of interneurons. Since the threshold value for firing of 
interneurons is rather small, we assume that when one pyramidal neuron fires, interneurons 
surrounding the firing pyramidal neuron immediately fire. Then, these firings of interneurons 
bring about inhibitory synaptic electric currents in all pyramidal neurons, because interneu- 
rons are connected to pyramidal neurons via inhibitory synapses. This inhibitory synaptic 
electric current /jp, which is independent of index /, is written as 



j k 



(6) 



where a-function 5jp(f ) is defined as 



V(0 = 



%,1 %,2 



for t < 0, 
for < f . 



(7) 



We set /iT = 1 /N for proper scaling. The constant Ajp is used to control the intensity of /jp, and 
the constants Tjp ^ and Tjp2 are always taken to be 10 [msec] and 1 [msec], respectively. The 
function 5'jp(f ) takes negative value so as to represent the inhibitory nature of the connection. 
Once note that we neglect the detailed dynamics of interneurons and simply assume that the 
inhibitory currents /jp are induced in all pyramidal neurons immediately after one pyramidal 
neurons fires [^7[]. 

The current 1^^^ ■ in /; is used to control initial firings of neurons. For the initial condition 
of the network, we set all state of neuron (v, , {w,-,}) to be at the stable fixed point of the 
dynamics (|l]) and with I- = 0. Since all neurons keep quiescent without any external 
stimuli, we use the pulsed form of the external electric current I^-^^j ■ to invoke initial firings, 
which act as a trigger for information processing of the present model. The detailed definition 
of /gj^j , will be given in section ^j. Note that the current /g^^y ■ is applied only in the beginning 
of the dynamics In the theoretical analysis below we always set I^^j^ ■ — because 

we focus on the stationary state in this analysis. 



2.2 Spike-timing-dependent learning rule 

We investigate associative memory neural network models that memorize multiple periodic 
spatio-temporal patterns of spike timing. P periodic spatio-temporal patterns to be memorized 
are generated randomly according to the equation 

.~f=^f+fcr, fc = ...,-2,-l, 0,1, 2..., ^ = 1,...,/^, A?, (8) 

with 

^f=^^r' (9) 

where Q is a. natural number controlling the degree of discreteness of spatio-temporal pat- 
terns, and random integer qf^ is chosen from the interval [0, Q) with equal probability. T 
denotes the period of the spatio-temporal patterns. We set T — 250 [msec] and 2 = 10 in 
what follows. 

Let us consider the problem of encoding the spatio-temporal patterns so that the net- 
works function as associative memory. The recent results of the electrophysiological experi- 
ments have revealed that the modification of a synaptic weight depends on the precise timing 
of presynaptic and postsynaptic spikes [ p^ ^ . Such modification of synaptic weight AJ 
is approximately written in the form 

Ay cx W(Af) 
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— (■ 



AAw.i_eAr/T^.2j forAf<0, 

(10) 

e-A'/%.i _ e-A'/%.2 \ for < Af , 



with 

= fpost-Ve' (11) 

where t^^^ and fpie denote spike timing of presynaptic and postsynaptic neurons, respectively. 
The asymmetric shape of the time window W(Af ) is described in Fig. |l], where we set t^^ ^ = 
25 [msec] and T^y 2 = 2.5 [msec] as we set in what follows. 

We encode the spatio-temporal patterns according to this spike-timing-dependent synap- 
tic plasticity. In our previous study [p^, we have introduced the learning rule 

where, to take account of the periodicity of the present spatio-temporal patterns s^, we define 
r-periodic function W{At) of the form 

W(Af) = £ W{At + kT) 

1 fe-A'/Vi _e-(^-A')/%,i _e-(r-A0/%,2" 



0<Af<r. (13) 

This learning rule is applied also to the present neural networks. As will be shown in the 
following sections, the spatio-temporal patterns encoded with this learning rule are retrieved 
successfully in the network dynamics (|l])-(^. 



3 Perfect retrieval state 

Here we investigate the stationary properties of retrieval state of the network in the limit of 
infinite number of neurons. In the present analysis, we focus on the retrieval state of the form 

tKk) = ^q\+kf, ^=...,-2,-1,0,1,2,..., /=l,...,Af, (14) 

where we suppose pattern 1 as the retrieved pattern. We term the retrieval state ( p4| ) perfect 
retrieval state since no spike timing is allowed to deviate from the encoded pattern in this 
retrieval state. Note that the period t in Eq. (|l4|) is different from the period T, which is 
assumed in generating the spatio-temporal patterns s^, that is, the period of the retrieval 
process is different from the period of the encoded pattern. In the present section, we aim 
to evaluate the period f, which determines the form of the periodic solution for the perfect 
retrieval state. The stability of the periodic solution is examined by a linear stability analysis 
in section Q 

One possible way to determine the period T is substituting Eq. ( p^ ) into Eqs. (^ and 
so as to obtain the periodic synaptic electric current = /pp ■ + /jp in the limit of ^ 00. 
Then, the current is evaluated as a function of f , and hence we obtain the periodic firing 
pattern of neurons as a function of T . Comparing the evaluated firing pattern with the 
substituted firing pattern ([l4|), we can determine the period T self-consistently ||l0|]. 

We follow the almost same scheme as above, although we make a slight detour for the 
convenience of the calculation below. We first consider the solution of the form 

ti{k)^t,{k), ^ = ...,-2,-1,0,1,2,..., /ef/,,, ^ = 0,...,e-l, (15) 
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where the set of indexes Uq is defined as 



(16) 



We term a cluster of neurons that belong to Uq sublattice q. In the solution (p^), neurons 
belonging to the same sublattice are assumed to behave in the same manner It will be shown 
that the dynamics has the solution of the form jl^ in the limit of ^ °o if f is 

finite [H]. We will evaluate the Q-boAy dynamics for these sublattices, which has important 
implication for the stability analysis in section Q After that, to this Q-body dynamics of 
sublattices, we substitute the solution of the form 

t*q{k)^^q + kt, A:=...,-2,-l,0,l,2,..., ^ = 0,...,e-l. (17) 

Then, we obtain the period f for the perfect retrieval state (p^. 

In the analysis below, we always assume finite P and finite Q. Asterisks are used to 
indicate the variables in the stationary state. 

3.1 Dynamics of sublattices 

In order to evaluate the dynamics of Q sublattices, we first evaluate the current /pp ■ in the 
limit of ^ oo under the condition (nsb. Assuming that neuron / belongs to sublattice q, we 
substitute Eqs. (|l2|) and (|l5]) into Eq. (Q). Then, we have 



'PP,i 



q'=ojeu^^ k 

— E— E ^ 



Ev f-vW 

k 



Ai>l q' JeU, 











E'^pp 

k 


t-t/k) 



Ipp 



EvE^p 



t-tq'ik) 



i e Uq, 



(18) 



where Nq denotes the number of members in Uq, and variables / , and W are defined as 



J ,=W 

m 



(P-l)W, 



(19) 



_ 1 



^=nE«'(g^) = n E ^[n^ 



Q 



(20) 



Equation ( p^ ) shows that the current /pp ^ is independent of index / in the limit of ^ oo. 
Thus, we define the sublattice variable /pp ^ as 



^pp 
Q 



E4/E5. 



pp 



t-tq,{k) 



ieUq, q^O,...,Q-l. 



(21) 



Following the same scheme, we rewrite the current /jp in Eq. in the form 

. A, 



^ip 



'IP 



EEv f-v(^) 



Q a' k 



(22) 
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Equations ( pT| ) and ( p2[ ) imply that the synaptic electric current Ij = Ipp ■ + /jp ■ depends 
only on q, that is, neurons belonging to the same sublattice obtain the same amount of synap- 
tic electric current. Therefore, the dynamics (|l])-(^ has the solution in which neurons be- 
longing to the same sublattice behave in the same manner, as we have assumed in Eq. (15). 
Such dynamics of sublattices is expressed as 



Vq = f[Vq,W 



with 



Iq — Jppg + V 



(23) 
(24) 

(25) 



where ^v^, 1*^^/1) represents the common state of neurons that belong to sublattice q. The 
common synaptic electric current in Eq. ( p5| ) is defined by Eqs. ( pl[ ) and (p^. 



3.2 Derivation of perfect retrieval state 

Let us find the periodic solution for the perfect retrieval state (|l7|) in the dynamics of sublat- 
tices (|3|)-(|5|). Substituting Eq. into Eq. @, we have 



fLfq,LS: 

^ q' k 



PP 



t - 



Q 



q +kf 



(26) 



where f -periodic function 5pp(f ) is defined as 

k 



1 



g '/"^ppj 



1 - g-^/'^pp.i 1 _ e-^/'^pp,2 
In the same manner, Eq. ( ^2| ) is rewritten as 

A, 



where f -periodic function 5jp(f ) is defined as 

k 

-1 



%,1 %,2 



l-g-^/V.i l_e-^Aip.2 



o<f < r. 



(27) 



(28) 



(29) 



Therefore, f -periodic solution for the perfect retrieval state {v*, {^9'}) dynamics 
of the form 



where 



I T* 

^'PP,f/+'lP- 



,«, ^=o,...,e-i, 



(30) 
(31) 

(32) 



As shown in Eqs. (^6|) and (|28|), /pp ^ and /j*p are functions of q and T, and also /* is a function 
of q and T. Hence, we can calculate the behavior of each sublattice as a function of q and T 
from the dynamics (pQ)-( 32 ). 
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Noting Eqs. @, and @, we obtain 

^^^^ 

It means that every synaptic electric current /* is identical, except that it exhibits the time 
shift according to q, and the behavior of all sublattices in Eqs. (^o|)-(^^ are evaluated from 
the time shift of sublattice 0. Therefore, we focus on the analysis of sublattice in what 
follows. 

We can calculate the behavior of sublattice in the dynamics (|30|)-(^2|) for the arbitrary 
value of T. In the stability analysis in section H, we will show that if the dynamics (Ub-P) 



has the stable perfect retrieval state, then the dynamics (|30|)-([32[) also has the stable periodic 
solution at f = T*, where T* denotes the solution of the period now under consideration. In 
almost every cases, for T that is sufficiently close to f *, sublattice in the dynamics (pO|)-(^) 
exhibits the periodic firing motion, and hence the spike timing of sublattice in the stationary 
state is written as [llOll 



tQ{k)=kT + r{T), ^=...,-2,-1,0,1,2... (34) 
On the other hand, from Eq. (p7|), we obtain the spike timing of sublattice as 

t^{k) = —0 + kt*=kr\ ^=...,-2,-1,0,1,2... (35) 

Comparing Eq. (^4|) with Eq. (^5|), we obtain the condition 

r(f*)=0. (36) 

As shown in the previous study [|ic|], we can easily evaluate the explicit form of the func- 
tion r{T) numerically by integrating the single body dynamics of sublattice in Eqs. (^0|)- 
(^^. Once we evaluate the explicit form of the function r{T), we obtain the solution T* from 
the condition (Bq). 



3.3 Optimal form of the time window W{At) to encode multiple spatio- 
temporal patterns 

In general, the properties of the network depend on the number of stored patterns P. We can 
encode a large number of patterns when the network exhibits the weak dependence on P. To 
see to what extent the properties of the network depend on P, we decompose /* = Ipp ^ + /j*p 
defined by Eqs. (E6) and (ES) into the form 



with 



' Q y 



MP 



(37) 



(38) 



Z*=App(P-l)W5pp(f) 



where 



1 



^ q'=Q 



pp 



Q 



(39) 



(40) 



We term Z* the crosstalk term since this term appears in Eq. (^7|) as a result of encoding mul- 
tiple spatio-temporal patterns. The function Spp{t) always takes the positive value. Hence, as 
P increases, Z* exhibits a increase or a decrease depending on the sign of W, until the perfect 
retrieval state breaks at the critical number of patterns P = P^. 

Let us take notice of W appearing in Eq. (^9|). The quantity W, which is defined by 
Eq. (20|), is the average of the time window W{At). For the time window W{At) defined by 
Eq. ( ly), one can easily show 

W = 0. (41) 
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In this case, the crosstalk term Z* vanishes, and hence we can encode the arbitrary number 
of patterns in the limit of ^ oo as far as P is finite. It turns out that the present form of the 
time window W{At), which is found in experiments, is of great advantage to reduce the size 
of W and also the crosstalk among encoded patterns. 



4 Stability of the perfect retrieval state 

Although we have derived the periodic solutions for the perfect retrieval state in the previ- 
ous section, it still remains unclear whether the derived periodic solutions are stable in the 
network dynamics (|l])-(H). In some cases, the derivation of periodic solutions in the previous 
section yields unstable solutions, and the network cannot settle into such unstable retrieval 
state. In the present section, we employ a linear stability analysis for the perfect retrieval 
state we have derived in the previous section. That is the application of Floquet theory, which 
yields an eigenvalue problem with the finite size of the matrix. 



4.1 Decomposition of the problem: stability of sublattices and stability 
of the perfect retrieval state in the dynamics of sublattices (PU[)-(P2[) 

In a linear stability analysis, infinitesimal perturbation is assumed in the initial condition, and 
then the time evolution of the deviation from the target solution is investigated to the first 
order in Taylor series expansion. When we apply Floquet theory to the present system, the 
spike timing of neuron / that belongs to sublattice q is written in the form 



t,{k)^t:,{k)+8t.ik), 



-2,-1,0,1,2..., ieUq, q^O,...,Q-l, 



(42) 



where we suppose pattern 1 as the retrieved pattern. We assume that the initial condition is 
correlated only with pattern 1 and the correlation with other patterns does not arise in the 
time evolution of the network dynamics, that is, we assume dtj{k) is correlated only with 



form 



1 , . . . , A^). Substituting Eqs. (|12|) and into Eq. (}^, we obtain /pp ■ (/ G of the 



Q / 



Q 



[q-q') 



^ q' k il>\ 




FF 




t;,{k)-5tj{k) 



^t*,{k)-5t.{k) 



(43) 



where we utilize the assumption that 5tj{k) is correlated only with ij = ^q^j [j = 1, . . . ,A^). 
Since Eq. ( p3| ) shows that /pp ^ depends only on q, we are allowed to define sublattice variable 



V,, = V,. = %^EVE7^ E s,,\t-t;,{k)-8t.{k) 



(44) 



Performing a truncated Taylor series expansion of Eq. (44), we have 



(45) 
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with 



c,' k 

where the derivative of 5'pp(r) is written as 



--t;,{k) dt^,{k), 



(46) 







for t < 0, 



-1 



e-'/^pp.2 I for 0<r, 



(47) 



"PP.l '■PP,2 \ "PP,! 



"PP,! 



and the sublattice variable 5tg{k) is defined as 



^Mk) = ^ L St,{k)- 



<? ieu„ 



Following the same scheme as /pp ^, we obtain the deviation of /jp in Eq. (^ as 



/jp — Iip + 5/jp, 



with 



^p = -^LLSip[t-t/k)]5t^,{k), 

^ a' k 



(48) 

(49) 
(50) 



where the derivative of 5'jp(f ) is written as 




for t < 0, 



S[p{t) = 



1 



1 



1 



e-'/V.2 I for 0<f. 



'IP,! '■1P,2 \ "-IP,! ''IP,2 

We represent deviation appearing in the state of neuron / by 



w,.. 



l = l,...,n, i eUg, q = Q,...,Q-l. 



(51) 



(52) 
(53) 



Noting Eq. (p4|), we safely replace I- — /pp ■ + /jp in Eq. (^ by sublattice variable Iq =tppg + 
/jp. Then, we perform a truncated Taylor series expansion of Eqs. (0)-(^ and obtain the 
dynamics of the form 



dv 



5v, + £ 



df 



dv 



Sw.,, + 81 (54) 



5w.,,, / = !,...,«,/ e f/^, ^ = 0,... ,g- 1 (55) 



with 



8lci = 5lpp^ + 5/ip, 



(56) 



where we introduce abbreviations such as 4^ 



dv 



From the definition of spike timing, we have v- \tq{k) + 5t^{k)\ = 6 {i & Uq), which yields 



5t,{k) = - 



where constant c is defined as 



(57) 



(58) 



Note that constant c is independent of q and k. Now we can evaluate the time evolution 
of (5vj, {5Wjy}) from Eqs. (p4|)-(^7|). To solve this dynamics we need to calculate 5lq = 
dipp ^ + 5/jp, in which Stq{k) are required at time t — t*{k) = . . . , — 2, — 1 , 0, 1 , 2 . . .) as 
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shown in Eqs. and (|^. We can evaluate 5ty{k) from {5v,- [t*{k)] } by use of Eqs. 
and @. 

It is a hopeless task to apply Floquet theory directly to the A^-body dynamics (p4|)-(p^) 
since that gives the eigenvalue problem with the infinite size of matrix. For the purpose of 
reducing the degree of freedom, we define the following sublattice variables 



K ^ ^" 



1 



in 



1=1,. ..,n, q = Q,...,Q-l. 



Then, from Eqs. ( p4| ) and 

5va = 



5w^, = 



, we have 



1 ieu„ 



dv 



df 



> OW; = 



dv 



If dw,, 



1 = !,...,«, q = Q,...,Q-l, 
where, from Eqs. @ and @, 5/^ in Eq. (|l]) IS written as 



(59) 
(60) 

(61) 

(62) 



t-t*Ak) 5t,{k) 



^IP 



g' k 

In addition, substituting Eq. (|7|) into Eq. (^8|), we have 

5t,{k) = llllJl 



II^Jp 



t-t^,{k) 51 /k). (63) 



(64) 



Now we obtain the g-body dynamics (^T])-(^^. Calculation of 5/^ requires 1 5t^, (A:) |, which 

are obtained from |5v^; t*, (k)] | t ogether with Eq. (|4|). To this Q-body dynamics we will 
apply Floquet theory in section ]4!2| 

The stability of the periodic solution in the Q-body dynamics (^T])-(^^ is the necessary 
condition for the stability of the retrieval state in the original dynamics^)-(|3]), but not the 
sufficient condition. Therefore, we must investigate the behavior of the following variables 



5v, 
5w., 



dv 



"1^ 



/ = !,...,«,/£ Ug, q = 0, 



(65) 

,2-1- (66) 



If the perfect retrieval state is stable, {dv-,\ dw;,]) converges into the fixed point (0, {0}). 
Subtracting Eqs. (^Tj) and (62 ) from Eqs. and ( p^ respectively, we obtain 



df 



5v, + I 



dv 



5v, + E 



I' 



/ = 0, . . . ,n, / G Uq, q = 0, 



(67) 

,G-1. (68) 



For the stable perfect retrieval state, the fixed point (0, {0}) is necessary to be stable in the 
dynamics and (Bsj). Note that deviations (5v,, {Sh",-,}) appearing in the dynam- 
ics Eqs. (|57|) and (|58[) do not interact with each other since this dynamics includes no in- 
teraction term like dly This stability problem is thus a single body problem, which is easily 
evaluated numerically. 

The stability problem of the perfect retrieval state in the dynamics (|lJ)-(§) is now decom- 
posed into two stability problems: the stability of the perfect retrieval state in the Q-body dy- 
namics (^T|)-(|64|) and the stability of the fixed point (0, {0}) in the single-body dynamics ( |67| ) 
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and (|68[). What are the imphcations of these two stabiUty problems? It is straightforward 
to see that the former problem is equivalent to the stability problem of the perfect retrieval 
state in the Q-body dynamics of sublattices (p3|)-(|2^. Hence, we conveniently call the for- 
mer problem the stability of the perfect retrieval state in the dynamics of sublattices. In the 
dynamics of sublattices (|23l)-(^5l) we neglect a distribution of spike timing of neurons in each 
sublattice, and this distribution of spike timing is treated in the latter problem. We thus term 
the latter problem the stability of sublattices. 

It is of interest that a truncated Taylor series expansion of Eqs. @ and with fixed /* 
gives the same stability problem as Eqs. ( |67| ) and (p8|). This result implies that if the periodic 

solution (^^'l*^^/}) stable in the dynamics (^0|)-(p2|), then the stability of sublattices 

are ensured. We evaluate the periodic solution ('^'9' 1*^^/}) numerical integration 

of the dynamics ( ^0| ) and (|3l| ), and hence it is impossible to obtain the unstable periodic 
solution of the dynamics ([3Cl|) and (|3l|). In other words, the numerically evaluated periodic 
solution ^v*, I'^J/I) is always stable, and also the stability of sublattices is always ensured. 
Therefore, further investigation on the stability of sublattice is unnecessary, and we focus on 
the stability of the perfect retrieval state in the dynamics of sublattices in the next section. 

4.2 Floquet theory for the perfect retrieval state in the dynamics of sub- 
lattices 

Here we apply Floquet theory to the Q-body dynamics (|6l|)-(|64}). In the evaluation of this 
dynamics, 5tq{k) are required at time t ^ t*{k) {k — . . . , —2,-1,0, 1,2 . . .). One may thus 
consider it convenient to define the vector 

5x,ik) - (5y,[f;(^)],5vv^^[f*(^)],...,5vv,„[f,*W]). (69) 

The vector 5x^(A;) represents the deviation at time t = t*{k). Since the vector 5xq{k) in- 
cludes the variable 5vg [tgik)] , we can calculate 5tq{k) from 5xg{k) by use of Eq. (|64|). Let 
us consider the problem of calculating the deviation 5xq(A:+ 1) from the past deviations 
8xq{k') {q = 0,...,Q-l, k' <k+l). 

The a-functions 5'pp(f ) and 5'jp(f ) give an infinite long-time influence after the activation, 
and the derivatives of these a-functions appearing in Eqs. ( |46| ) and also have an infinite 
long-time influence. It means that long past deviations 5tq{k') and also 5xq{k') are necessary 
in the evaluation of the present value of 5Iq. It is again a hopeless task to consider Floquet 
theory based on the vector 5xq{k) since that still gives an eigenvalue problem with the infinite 
size of matrix. 

For the further reduction of the size of matrix, we define the variables 



A 



A, 



-\t-e,{k')-st^,{k')]ix^^,_ 



h^. - ^LV L --^-^-^^^ (71) 



1 



^ q' t'i{k')<t ^IP,1 ^IP,2 



1 



-[(-r;,(A-')-5r^,(i')]/V,2 



'4,, = E " _ • - oy, 



1 



Then, for the specific form of a-function (Q), we can rewrite Eq. (^6|) as 
5^pp.o = -^LJoq'Spp['-';'ik)]St^'ik) 
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51^ [fo* {k)]e- i'-'o ('^l] /"pp.1 + S/j^oK {k)]e- i'-'o ('^^l /"pp,2 , 

t^{k)<t<t^{k+l), (74) 



where 



m)-'*Jk')\/rpp,, 



5/:,oKW] = -^EV L ^ ^^Ifik'), (75) 

';,(':')<'o*W '^PP,l(^PP,l-'^PP,2j 

5/2,o[^o(^)] = -^LV L 7 Y^VC^')- (76) 

9' ';,('^')<'o*W%P,2(%P,l-%P,2j 

In the same way, we rewrite Eq. ( pO[ ) as 

fo*W<f<fo(fc+l), (77) 



where 



4 .-['oW-'J/C^Ol/Tip,, 
5/3oK(^)] = --^L L 7 -St/k'), (78) 



',1 '■IP,2 



In Eqs. (^ and (0), 5/pp q and 5/jp (fg (A:) < f < fg 1 ) ) are evaluated only from 1 5?^, (A:) | 

and 1 5/., q [f^ (A)] | . Solving the dynamics ( ^ and ( |62| ) with S/q = 5/pp q + 5/jp defined by 
Eqs. ( [74I ) and ( [tt] ) under the condition (5v()[fQ(A)], {5wQ^,[fQ(A)]}), we obtain the next devia- 
tion {5vq [t^{k +!)],{ 5wqi, [f^ + 1 )] } ) as a function of 1 5t^, (k) | , 5vq [t^ (k)], { 5 Wg;, [t^ {k)]}, 

and 1 5/^/ o[fo(^)]|- Hence, we are allowed to define the functions 

5voK(fc+l)] = /?({57,,(fc)},5vo[foW],{5vvo,,[fo*W]},{5/,,oKW]}), (80) 

[fo* (fe + 1 )] = 5, ( { 5F^, (k) } , 5 vo K m,{ Swoi' K (^)] } , { 5/.,o [fo (k)]}), 

1 = 1,. ..,n. (81) 

Because of the form of the dynamics ( |6l| ) and (^), we obtain 
^({57^,},5v,{5w,,},{5/,}) 

S,[{dt^,},5v,{5w,,},{dl,}) 

55, 55, „ dS, „ 55, 



/=l,...,n. (83) 

Note that every coefficient in Eqs. (^) and (|3]) is a constant, which is independent of 
({57^,},5v,{5vv,},{5/,}). 
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Meanwhile, from Eq. (^5]), we obtain 

5/i,o[fo(^+l)] 



A .-['o('^+i)-OW]/Tpp.i 



(84) 

It means that 5/j Q[fQ(^+ 1)] is a function of |5f^,(A:)| and 5/j o[fo(^)]- We obtain the sim- 
ilar relation for the rest of 1 5/,, ^ [fg + 1 )] | , and they are also functions of | St^^, (k) | and 

Now, we define the vector 

5x, {k) = (5y, [t; (k)] , [t; (k)],..., dw,„ [t^ (k)] , dl,^ [f * {k)],..., dl,^ [t^ {k)] ^ 

;t=...,-2,-l,0,l,2..., ^ = 0,...,e-l. (85) 

Then, noting Eq. (^4|), we can summarize Eqs. (|o|)-(|^, and so on in the form 

e-i 



5xo(^+l)= 5xc,{k)+B5xQ{k), 



(86) 



where the definitions of the matrices and B are given in Appendix ^. Furthermore, we 
define the vectors 



(87) 
(88) 



5X(0) = (^5xg_i(^),5xg_2W,---,5xoW 
5X(1) = (5xo(A:+l),5xg_i(^),...,5xi(^) 

Then, the relation between 5X(0) and 5X(1) is written as 

5X(1)=M5X(0), 

where 



(89) 



M: 



/ Agj Ag_ 
E 
E 



Ai B \ 





(90) 



V E / 

Because of the symmetrical properties of the present system, we have 

5X(2) = M 5X(1) = 5X(0), 

where 

5X(2) = (Sxi + 1 ) , 5xo(^ + 1 ) , 5xg_ i (^) , . . . , 5% {k) 
Following the same scheme, we obtain the vectors of the form 



5X{n)^M" 5X(0), 



1 <n, 



(91) 
(92) 

(93) 



where 5X(n) represents the deviations in the future. 

Now the stability problem of the periodic solution is reduced into the eigenvalue problem 
with the finite size of the matrix M. As will be shown in Appendices ^ andO, we can easily 
evaluate the matrix M numerically since the coefficients i n Eq s. ( ^2| ) and (^Sfare obtained by 
numerical integration of the single-body dynamics ( IH 15[ ). In general, the matrix derived 
in Floquet theory always has the eigenvalue Aj = 1 with the eigenvector corresponding to the 
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time shift in the periodic solution. The matrix M also has the eigenvalue Aj = 1 with the 
eigenvector 

5Xo = (5xo,...,5xo) (94) 

with 

5x, = (v5 [t*o (k)] , vvSi K m , . . . , vv5„ K (k)] , n.o [tmi ■ ■ ■ , 4*.o k m) , (95) 

where we define |/*o[foW]| by substituting dt^^,{k') = (/t' = ...,-2,-1,0,1,2..., q' = 

0, . . . ,Q — I) into the derivatives of Eqs. ([70[)-([73|). If the periodic solution is stable, the abso- 
lute value of other eigenvalues \X,„\ ( 1 < m) must be less than 1 . Therefore, we can determine 
the stability of the perfect retrieval state by numerical computation of the eigenvalues of M. 

In the following sections, we will apply the present analysis to evaluate the stable perfect 
retrieval state for the various value of parameters. As will be shown, the present analysis is 
powerful enough to draw the phase diagrams. 



5 Retrieval process 

In this section, we illustrate the typical behavior of network in the process of memory re- 
trieval. In what follows, we always assume Q = 10, which brings about discrete type of firing 
pattern of memory retrieval. For the initial condition of the network, we set all state of neuron 
(vj, {w,/}) to be at the stable fixed point of the dynamics (|l]) and (H) with ~ 0. To evoke 
the retrieval of pattern 1 , we give the external stimuli of the form 

l^^^ . = I ^EXT^ (f - %^?/) for < f < AextText' (95) 
[ otherwise. 

where 5(f) represents the delta function, and the parameters Agj^.^, 7"^^,^, and flgxx chosen 
so that the initial part of the pattern 1 is forced to be retrieved. In the present study, we set 
"^EXT ~ 3*^' ^EXT ^ ^^'^ ^Exr < 0.1. Note that the external stimuli I^^j ■ is applied only in 
the beginning of the network dynamics. 

In Fig. §(a), we describe the result of the numerical simulation with P = 3 and = 8000. 
The initial firings of the neurons are evoked by the external electric current /g^^ j, while 
other firings are brought about by the synaptic electric current /pp ■ +/ip. The firing pattern 
in Fig. ||(a) , which looks like vertical bars, indicates the synchronized firing of a numerous 
number of neurons. Since it is difficult to see whether the retrieval of pattern 1 is realized in 
Fig. ||(a), we replot the same result of numerical simulation in Fig. §(b), where the vertical 
axes is set to represent . In this figure, we clearly see the successful retrieval of pattern 1, 
in which the neurons belonging to the same sublattice exhibit synchronized firing. 

The dynamical behavior of the neuron with qj = is described in Fig. ||(a) . After the 
transient behavior, the neuron settles into the stationary state, where the neuron exhibits pe- 
riodic firing. In Fig. ||(b), we describe the periodic solution for retrieval state obtained from 
Eq. (^). In order to examine the stability of this solution, we calculate the explicit value of 
the matrix M numerically. In the present case, the largest absolute eigenvalue is 1, and the 
theoretically evaluated perfect retrieval state in Fig. ||(b) is stable. The good agreement be- 
tween Fig. ||(a) and (b) implies the validity of the present analysis. It is also worth noting that 
the theoretical result in Fig. p|(b) is independent of P because of Eq. (p4|). We set P 3 in the 
numerical simulation in Fig.^a), and this result of numerical simulation is well explained by 
the P-independent solution in Fig. ^b). We will see the same result of numerical simulation 
even with the larger value of P, as far as P/N is sufficiently small. 

6 Phase transition due to the loss of the stability of the per- 
fect retrieval state 

Here we investigate the effect of inhibitory synaptic electric current /jp, which is controlled 
by Ajp. From Eq. (36), we obtain f /2 as a function of Ajp, which is plotted in Fig. ^a) . As 
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Ajp increases, the period of the retrieval process T becomes longer since each neuron obtains 
a large amount of inhibitory synaptic electric /jp with the large value of Ajp. Figure ^b) 
describes the absolute eigenvalues of the matrix M. Size of the matrix M is 80 x 80, and the 
largest two absolute eigenvalues are plotted in Fig. ^b) . With Ajp^cSOO the largest absolute 
eigenvalue is 1, while it exceeds 1 with 500;£Ajp, that is to say, the stability of the perfect 
retrieval state is lost beyond the critical point Ajp ~ 500. 

To observe this phase transition in numerical simulations, we calculate the inter spike 
intervals (ISIs) of all neurons changing the value of Ajp, as described in Fig. ^c). Note that 
the ISIs we calculate here is based on spike timing of all neurons. When neuron / and neuron j 
fire sequentially at time and respectively, we calculate the time difference f ^ — to obtain 
the ISIs of all neurons. By means of these ISIs, we can evaluate the gaps in spike timing 
appearing in Fig. ^a), which corresponds toT /Q. The theoretical result in Fig. ^a) and (b) 
explains ISIs in Fig. ^(c) well, although we see some fluctuations due to the finite number 
of neurons near the critical point Ajp. In Fig. ^d), we calculate the ISIs in the dynamics of 
sublattices (^3])-(|5|), in which we have taken the limit of the infinite number of neurons. The 
theoretically evaluated critical point A[p explains the loss of the stability observed in Fig. ^d) 
with a high degree of precision. 

In Fig. H we draw Ajp-App phase diagram, which is evaluated by the theoretical analysis. 
We find the stable perfect retrieval state in the region represented by PR. As App decreases, 
the range of Ajp for the stable perfect retrieval state becomes narrower since a large amount 
of /pp is required for the successful memory retrieval under the strong inhibition. 



7 Two separated perfect retrieval phases appearing with 
the slow a-f unction S^^{t) 

In the previous sections, we assume 5'pp(f) with Tpp j — 3 [msec] and Tppj — 0.3 [msec] as 
well as 5'jp(f) with Tjp j — 10 [msec] and Tjp2 ~ ^ [msec], where 5'pp(r) decays much faster 
than 5jp(f). In order to examine the role of the decay time constants in a-functions, we inves- 
tigate the case of the slow a-function 5pp(f ) with Tpp j — 20 [msec] and 5pp 2=2 [msec]. For 
this slow a-function 5'pp(r), we describe Ajp — App phase diagram in Fig. ^. The distinctive 
feature of this phase diagram is the perfect retrieval phase appearing in the region with the 
large value of Ajp. In the case of the fast a-function 5'pp(f ), the strong inhibition with the 
large value of Ajp tends to suppress the perfect retrieval, as described in Fig. ^j. Nevertheless, 
in Fig. ^, we see two separated perfect retrieval phases in the region with 40000 ;<App;; 70000, 
while these two retrieval phases merge with each other in the region with 70000<^App. 

One example of the retrieval process in the region with the large value of Ajp is illustrated 
in Fig. 0. As a result of the large value of Ajp, the neuron obtains a large amount of the in- 
hibitory electric current /jp, which oscillates with the period f /Q. In the retrievalprocess, Q 
sublattices emerge exhibiting synchronized firing of neurons, as described in Fig. |(a). When 
one firing of sublattice occurs, all neurons obtain a large amount of the inhibitory synaptic 
electric current /jp. Then, neurons in the next sublattice cannot fire until this inhibitory elec- 
tric current decays with the time constant Tjp j . In this way, the oscillatory inhibitory electric 
current /jp regulates the spike intervals of sublattices, and hence the memory retrieval with 
the long period T is realized. The long-time influence of the slow a-function 5pp(f) is indis- 
pensable for this memory retrieval since the time gaps of firings of sublattices (i.e., t /Q) are 
considerably large. 

In Fig. ^ (a) and (b), we describe the result of our analysis with App = 65000, where we 
see two separated perfect retrieval phases. Near the boundary of the retrieval phases, Eq. ( p6| ) 
yields two different perfect retrieval states, which are indicated by 's' and 'u' in Fig. ^ (a). 
As described in Fig. ^ (b), the largest absolute eigenvalue of the matrix M for the state u 
exceeds 1, while that for the state s always takes 1. This result implies that the state s is stable 
and the state u is unstable. The ISIs observed in the numerical simulations are plotted as a 
function of Ajp in Fig. ^ (c). To obtain these ISIs, we slowly change the value of Ajp both 
from Ajp — and from Ajp — 2000. In the present case, neurons cease firing between the 
critical points Ajp(l) and A^(2). The ISIs observed in the dynamics of sublattices (^3|)-(^) 
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are also plotted in Fig. ^ (d). The phase transitions observed in Fig. ^ (c) and (d) are well 
explained by the theoretical analysis in Fig. |](a) and (b). 

In order to investigate more details about decay time constants, we describe Ajp — Tpp j 
phase diagram in Fig. ^, in which we fix Tpp, = 0.1 Tpp j. In the region with the long Tpp , 
we find two separated retrieval phases. In this case, stable and unstable solutions of Eq. ( P6| ) 
are found only inside the perfect retrieval phase, as described in Fig. |[ One might thus 
conceive that the stability analysis is not necessary for the purpose of determining the phase 
boundary. However, with the short Tpp j , we find the unstable solution of Eq. ( |3^ outside the 
perfect retrieval phase, as described in Fig. ^ The stability analysis is hence indispensable to 
determine the boundary of the perfect retrieval phase, particularly with the short Tpp ^ . 



8 Discussion 

We have investigated associative memory neural networks of spiking neurons memorizing 
periodic spatio-temporal patterns of spike timing. In encoding the multiple spatio-temporal 
patterns, we assume the spike-timing-dependent synaptic plasticity with the asymmetric time 
window W{At) in Fig. |l[ Encoded periodic spatio-temporal patterns of spike timing are 
reproduced successfully in the periodic firing pattern of neurons in the process of memory 
retrieval. In this retrieval process, Q sublattices (clusters of neurons) exhibit synchronized 
firing, and the oscillatory inhibitory electric current /jp, which is supposed to come from 
interneurons, regulates the spike timing of sublattices. 

In order to investigate the stationary properties of the system, we have derived the periodic 
solution for the retrieval state analytically in the limit of infinite number of neurons. From 
this analysis, we have shown that if the average of the time window W(Af) takes the value 
of zero, the crosstalk among encoded patterns vanishes. This result implies that the present 
form of the time window W{At), which is found in experiments, has the great advantage in 
encoding a large number of spatio-temporal patterns. 

To elucidate the stability of the derived periodic solution we have employed a linear sta- 
biUty analysis. In this linear stability analysis we have to evaluate the time evolution of 
infinitesimal deviation so as to obtain the matrix for Floquet theory, although the naive appli- 
cation of Floquet theory yields infinite size of matrix. In order to reduce the size of matrix, we 
have employed some decomposition of the stability problem, by which the original stability 
problem with neurons is reduced into the stability problem with Q sublattices. Then, to take 
account of the infinite long-time influence of a-functions, we have introduced the variables 
, which enable us to obtain the finite size of matrix M for Floquet theory. The explicit 



form of M is computed by the numerical integration of the single-body dynamics (111)-(115), 
and the stability of solutions are evaluated from the eigenvalues of M. 

Based on these method of analysis, we have investigated the stationary properties of 
retrieval state in the case of the fast a-function 5pp(f) with Tpp ^ = 3 [msec] and Tpp, — 
0.3 [msec]. In Fig. ^a), we have obtained periodic solutions for the retrieval state for various 
value of App by solving Eq. (^). Then, in Fig. ^b), we have employed the stability analysis 
of these periodic solution to obtain the critical point Ajp . The phase transition observed in 
the numerical simulations in Fig. ^c) and (d) is well explained by this critical point Ajp. The 
condition for the successful memory retrieval is summarized as Ajp — App phase diagram in 
Fig. I 

Meanwhile, with the slow a-function Spp{t) with Tpp j ~ 20 [msec] and Tpp 2 = 2 [msec], 
we have found two separated retrieval phases, as shown in Fig. ^. The behavior of neurons 
in the memory retrieval with the large value of Ajp is described in Fig. ^ where we see the 
large size of oscillatory inhibitory synaptic electric current /jp regulating the spike timing of 
neurons. The result of the theoretical analysis is illustrated in Eq. (Q), where the stability 
analysis is used to chose the stable solution from the multiple solutions of Eq. (^6|). 

The heart of the present stability analysis lies in the exact reduction of the size of the 
matrix for Floquet theory. Since Q sublattices arise in the stationary state, we have to evaluate 
the matrix with the dimension of (1 +n + 4-)Q, where 1 + n corresponds to the degree of 
freedom of the neuron dynamics /(vj,w,[, . . . and gi(^V',Wji, ... ,Wj„) (/ = 1, ...,«), and 
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additional 4 is required to evaluate the infinite long-time influence of a-functions Spp{t) and 
5jp(f). In the present study we set Q = 10 and n — 3, which yields the matrix with the 
dimension of 80. Although one might conceive that the size of this matrix is somewhat large, 
the critical points obtained from this matrix well explain the result of numerical simulations 
with a high degree of precision, as demonstrated in Figs. ^ and||. In other network models 
[^^Kl, only a few number of sublattices emerges and the size of matrix becomes small. 

The spatio-temporal patterns to be memorized are assumed to be periodic in the present 
study for ease of analysis. It is worth noting that learning rule based on the spike-timing- 
dependent synaptic plasticity is applicable to a wide class of spatio-temporal patterns of spike 
timing. The periodicity of spatio-temporal patterns is not crucial, and it is almost obvious that 
spike trains generated by independent Poisson process are also well encoded by use of the 
time window W(Af). In the case of Poisson process, the firing rate in Poisson process must be 
adequately low since the refractoriness of neurons is expected to prevent retrieval of Poisson 
trains with high firing rate. 

It is of interest to consider the effect of noise in the present model. With the large value of 
Ajp, neurons obtain the large size of oscillatory inhibitory electric current Ipp as described in 
Fig. 0, and the effect like stochastic resonance is expected to occur in the presence of noise. 
The evaluation of the effect of noise, however, seems to be difficult in the present scheme of 
analysis since we are required to calculate the distribution of spike timing of neurons in this 
evaluation. 

With the large Ajp and the short Tpp j the basin of attractors for spatio-temporal patterns 
are found to be narrow in the present model (data not shown). In the initial condition, the 
inhibitory synaptic electric current /jp is taken to be the value of zero. With the short Tpp j , 
firing of the first sublattice, which is induced by I^y^j brings about firing of the second 
sublattice immediately since some accumulation of inhibitory synaptic electric currents /jp 
are necessary to control the next firing. For these reasons, a first few firings of sublattices 
take place quite rapidly. These rapid firings of sublattices give rise to too much accumulation 
of inhibitory electric current /jp, and then terminate firings of all neurons. The core of the 
problem in this phenomenon is too rapid firings of interneurons. To avoid this problem, 
more sophisticated modeling of interneurons is needed so as to realize adequate control of 
interspike intervals of interneurons. When we assume that interneurons exhibit periodic firing 
independently of pyramidal neurons, the inhibitory synaptic electric current /jp takes the form 

Ap=VLV(f-^71p), (97) 

k 

where Tjp represents the period of firings of interneurons. We can investigate the case of 
this periodic inhibitory electric current I^p following the almost same scheme of the present 
analysis. 

Finally, we discuss the biological implication of the present study. The result of the 
present study strongly suggests the possibility of the concept of temporal coding, in which 
information is assumed to be processed based on spike timing of neurons. The question then 
arises about where we can find these kind of information processing in the real nervous sys- 
tem. It is well known that the hippocampus is the important tissue for short term memory. In 
CAS region of hippocampus, we see dense recurrent connections among pyramidal neurons, 
and hence the short term memory is thought to be stored in the CA3 region of hippocampus. 
Memory stored in hippocampus should be transfered into other regions such as neocortex 
so that it is stored as the long term memory. Recently, some experimental results begin to 
suggest that this memory transfer process takes place when sharp waves (SPW) appear in 
hippocampus [^8|]. In SPW, fast periodic firings of interneurons (^^200 [Hz]) bring about 
oscillatory inhibitory synaptic electric currents in pyramidal neurons, and these oscillatory 
electric currents regulate occasional firings of pyramidal neurons [|o|]. Nadasdy et al. have 
investigated these occasional firings of pyramidal neurons and revealed that repeating firing 
patterns of pyramidal neurons are present in SPW [p^. These results of experiments indicate 
that spike timing of pyramidal neurons of SPW represent some kind of memory that should be 
transfered into neocortex. Also in the gamma oscillation, we observe oscillatory inhibitory 
synaptic electric currents due to periodic firing of interneurons, although its frequency is 
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somewhat low (20-80 [Hz]). Buzsaki et al. have hypothesized that the firing patterns of pyra- 
midal neurons in the gamma oscillation are stored in the recurrent connections of the CAS 
region of hippocampus, and then these stored firing patterns are replayed in the firing pat- 
terns in SPW in the time compressed manner Our theoretical model explains these time 
compressed replay of firing patterns. 

Some aspects of our theoretical model are, however, still biologically implausible. For 
example, the learning rule ( [l2| ) gives either negative or positive synaptic synaptic weights 
by chance although synaptic weight among pyramidal neurons are found to be positive in 
experiments. More precise modeling of interneurons might be needed to acquire a deeper 
understanding of the time compressed replay of firing patterns. Solving these problems will 
be part of our future study. 

A The Hodgkin-Huxley equations 

The Hodgkin-Huxley equations are the ordinary differential equations with four degrees of 
freedom, which have been developed to describe the spike generation of the squid's giant 
axon 1^]. In the present study, for the dynamics /(v,Wj, ... ,w„) and (v, Wj, . . . , w„) (/ = 
1 ,...,«), we assume the Hodgkin-Huxley equations of the form 

Cm / 

with 

«2 

^2 
«3 

where v represents the membrane potential, and and the activation and inactivation 
variables of the sodium current, and the activation variable of the potassium current. 
The values of parameters are — 50 [mV], Vj^ — —11 [mV], — —54.4 [mV], = 
120 [mS/cm^], = 36 [mS/cm^], = 0.3 [mS/cm^], and C,„ = 1 [^uF/cm^]. 



•••'^*'3) = ^A'«^*'2W'l(^Af„-^')+5/f^3(^'jf-^) 



+?l(v'l-v), (98) 

...,W3) = ai(l-wi)-j3iw,, (99) 

...jWg) a2(l - wj) -jSjWj, (100) 

...jWg) = a3(l - W3) -jSgWj (101) 

= 0.01(10- v)y|exp(^l^j-l|, (102) 

= 0.125exp(-v/80), (103) 

= 0.1(25-v)/|exp(^^^-l|, (104) 

= 4exp(-v/18), (105) 

= 0.07exp(-v/20), (106) 

. l/fexpf^LjlVll. (107, 
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B Definition of tlie matrices and B 

From Eqs. (||),(|85)-(0, and so on, in Eq 



1 dR 



1 ds„ 



c 






(f-^f/ej/Tpp J 












■"^PP.l^^PPi) 


Ajpg 


-if/Q)/^-iPi 


cgTip 1 ( 




-( 

-A,pe 


r-,f/e)/T,p2 



V '•'2v.2(v,i V.2 
In the same way, we obtain B in Eq 



|8^)as 



B 













dR 



/ I dR . dR dR 
1 aS, dS, dS, 



i_dSjj_ I dS„ 

"e d{Stg) d{Sv) 



dS„ 



cQx, 



PP.H'PP,! ^,2) 



-r/T, 



PP,2 



e2i^PP.2('''pp.i '^vv.i) 



cQt, 



IP.l I 'IP.l 



with 



\ '^■2V.2(v.i V,2 










IS written as 

\ 











dR 
d{Sw„) 
dS, 



q^l,...,Q-l. 



(108) 



d{Sv) dl^ 5{5W„) 



d(Sw„) 








(109) 



c = 



dR 


dR 


dR 


dR 


\ 


d{Sl,) 
dSi 


d{Sl,) 
dSj 


a (5/3) 

dSj 


d(Sl,) 


a (5/,) 


d{Sl,) 


a (5/3) 


d{Sl,) 




dS„ 


dS„ 




dS„ 




d{Sl,) 


d{Sl,) 




a (5/3) 



d{Sl,) 











































J 



(110) 



V 

In our analysis, we have to evaluate the eigenvalues of M numerically, and hence the 



numerical evaluation of Ag and B is required. The coefficients appearing in Eqs. ( |82| ) and (g3| 
are evaluated in Appendix ^ Except for b^^, . . . ,b^^, elements in Ag and B are determined 
by use of the values of coefficients obtained in Appendix ^, while we set , . . . , so that 
M has the eigenvalue Aj = 1 with eigenvector (|94|). 
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C Numerical evaluation of the coefficients in the functions 

R({5t,},dV,{dW,,},{BI,})andS,{{St,},5v,{5w,},{5I,}) 
{/ = !,...,«) 

In order to evaluate the coefficients in the functions R[...] and . .] (/ = 1, . . . ,n), we con- 
sider the single-body dynamics of the form 



/ = l,...n 



with 



where 



'pp,o 



and 



^0 ^ ^PP,0 +Ap' 



(111) 

(112) 
(113) 



4pp V f S 



t t;,{k) 5t^,{k)\ + (/r,oK«] + 5A_oK(^)])^"['"''*''']/''" 



/l.oKW] + 5/2nKW])^-['-''*f"]/'P''-^ m < t < t*o{k+l) (114) 



t-t;,{k)-5t^,{k) 



+ (/4*,o[fo(^)] + 5/4,oKW])^^^'"''**"]/''^^ m < t < t^{k+l). (115) 
From Eqs. (^-(^, we obtain the explicit form of |/*o[f() (^)] | as follows 

g-(f-./f/e)/Tpp, 



A*oK(^)] 



'^PP.l %'P,2 I ' 1 



)(l_e-?App, 
-(f-,'f/e)/Tpp_j 



Q Y'"' 

A 

2pp^j ^ 

Q q' (Tpp,l-Tpp,2)(l-. 

Ajp^ -g-(^-'^'y/e)/v.. 

Q ci' (%,i-%,2)(l-^ 
Ajp ^ ,-{T-c,'TlQ)lz 



2" %,l-%,2l'l-^-'/^- 



We solve the dynamics ( 1 1 1 )-( 1 15 ) under the condition 



"01 



I ~ 1, . . . ,n. 



(116) 



(117) 



(118) 



(119) 



(120) 
(121) 



Then, we obtain the functions 



1=1,. ..,n 

It is straightforward to show 

R{{dt,},5v,{dW,,},{dI,}) 



(123) 
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^^ F{{£dl,},£dv,{£dw,},{£dl,})-Vl[t*,ik+1)] 
e^O £ ' 

5,({5F,},5v,{5w,,},{5/,}) 

^.^ G,{{£5t,},£5v,{£5w,,},{£5l,})-w*,i[t*,{k+l)] ;^i_^„(i25) 
£^0 e ' 1 • • • 1 J 

where < £5Jq . The expHcit value of F[. ..] and G/ [. . .] is easily computed by the numerical 



integration of the dynamics (1 1 1)-(1 15). We obtain 7? [. . .] and S'J. . .] by evaluating F[...] and 
G;[. . .] for sufficiently small £. 

Once we obtain R[. . .] and . .], we can easily evaluate the coefficients appearing in 
Eqs. (||) and (|3|). For example, substituting ({0}, 1 , {0}, {0}) into Eq. we have 

/?({0},1,{0},{0})^^. (126) 

Hence, is calculated from /?({0}, 1, {0}, {0}), which is computed by Eq. ( |124 ) with 
sufficiently small £. In the same manner, we obtain every coefficient in Eqs. ( p2| ) and (83|). 
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Figure 1: The shape of the time window W{At) with ^ = 25 [msec] and % 2 = 2.5 [msec]. 
The modification of synaptic weight is written as A/ oc W{At), where At = t^^^f — tpie denotes 
the time difference between the postsynaptic and presynaptic firings. 
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Figure 2: The retrieval process of pattern 1 observed in the numerical simulation with App = 
30000, Tpp J = 3 [msec], Tpp j = 0.3 [msec], Ajp = 250, Tpp ^ = 10 [msec], Tpp 2 = 1 [msec], 
P = 3, and = 8000. (a) Spike timing of neurons are plotted by closed circles as a function 
of time. The initial firings of the neurons are evoked by the external electric current ^ext!' 
while other firings are brought about by the synaptic electric current Ipp^+I^p. (b) Setting 
the vertical axis representing we replot the same result of the numerical simulation. The 
neurons belonging to the same sublattice exhibit synchronized firing, which are expressed as 
overlapped closed circles. 
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Figure 3: (a) The behavior of a neuron with qj = observed in the numerical simulation in 
Fig. |[ The membrane potential v. and the synaptic electric current I- are plotted as a function 
of time, (b) The result of the theoretical analysis for the stationary state of the neuron, which 
explains the stationary behavior of the neuron in (a). 
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Figure 4: (a) f /Q obtained from Eq. ( |36| ) is plotted as a function of Ajp. (b) The largest two 
absolute eigenvalues of the matrix M are plotted as a function of Ajp. One of the absolute 
eigenvalues exceeds 1 beyond the critical point App ^ 500, which is represented by the verti- 
cal lines in all four figures. The perfect retrieval state we have evaluated in (a) is stable only 
below the critical point App. (c) Changing the value of Ajp we observe the ISIs of all neurons 
in the numerical simulations with P —\ and = 8000. See text for the definition of the ISIs 
we calculate here, (d) We observe the ISIs also in the dynamics of sublattices (p3|)-(p5|). The 
critical point Ajp evaluated in (a) and (b) explains the phase transition observed in (c) and (d) 
well, although we see some fluctuations near Ajp owing to the finite number of neurons in the 
case of (c). Beyond the critical point Ajj, the network settles into another stationary state. The 



value of parameters are Tpp j = 3 [msec], Tpp 2 = 0.3 [msec], App = 30000, Tj 



and T, 



IP,2 



1 [msec]. 



10 [msec]. 
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Figure 5; Ajp — App phase diagram obtained by the theoretical analysis. The stable perfect 
retrieval state is found in the region denoted by PR. The value of parameters are Tpp j ~ 
3 [msec], Tppj — 0.3 [msec], Tjp j — 10 [msec], and Tpp j = 1 [msec]. 
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Figure 6: The same as Fig. H, except that Tpp j — 20 [msec] and Tpp 2=2 [msec]. 



29 



Masahiko Yoshioka 



(a) 




200 400 600 800 1000 
t [msec] 



(b) 



> 
> 

C 
CO 

E 
o 

i 



100 






1 








V 


50 





























-50 


























100 











200 400 600 800 1000 
t [msec] 



Figure 7: The retrieval process of pattern 1 in the case of the slow a-function 5pp(f) with 
Tpp J = 20 [msec] and Tpp 2=2 [msec] under the strong inhibition with Ajp = 1500. (a) Spike 
timing of neurons observed in the numerical simulations with P — 3 and A' = 8000 are plotted 
as a function of time. Note that the vertical axis represents qj . (b) The membrane potential v, 
and the synaptic electric current Ij are plotted as a function of time for a neuron with qj = 0. 
The value of parameters are App = 65000, % j = 10 [msec], and % 2 ~ ^ [msec]. 
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Figure 8: (a) For the case of the slow a-function 5pp(f ) with Tpp j = 20 [msec] and Tpp 2 = 
2 [msec], we plot T jQ obtained from Eq. ( |3^ ) as a function of Ajp. Near the boundary, we 
find two solutions of Eq. (^), which are represented by 's' and 'u'. (b) The largest absolute 
eigenvalue of the matrix M is plotted as a function Ajp. The eigenvalues for state s and state u 
are represented by Aj and A„, respectively. This result implies that state s is stable, while 
state u is unstable. The two critical point Ajp(l) and Ajp(2) obtained from the theoretical 
analysis (a) and (b) are represented by the vertical lines in all four figures, (c) The ISIs 
observed in the numerical simulations with P — \ and = 8000. We slowly change the 
value of Ajp both from Ajp — and Ajp = 2000. (d) The ISIs observed in the dynamics of 
sublattices (p3|)-(p5|). The phase transitions observed in the numerical simulations (c) and 
(d) are well explained by the critical points Ajp(l) and AJp(2). In (c) and (d), all neurons 
cease firing between Ajp(l) and Ajp(2). The value of parameters are App = 65000, Tjp ^ = 
10 [msec], and % 2 ~ ^ [msec]. 
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Figure 9: Ajp — Tpp , phase diagram obtained by the theoretical analysis, where we fix Tpp ^ = 
0. 1 Tpp J . The value of parameters are App = 65000, % j = 10 [msec], and % 2 = ^ [msec]. 
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